clear all
version 16
cd "_______________"
set more off

// Stata 16 has improved in optimization, at least with our empirical example. 
// Most city years converge when Stata automatic starting values are used.
// A few cities do not converge with Stata automatic starting value. In most cases changing starting value to estimates from the next converged years of the same city will do the trick. 
// There are some difficult city/years. 
// id 73 year 2001 converged with starting value 0.34, -1.86, 2.26, 1.98.
// id 45 year 2002, 2003 converged with starting value 2.78, -1.73, -1.14, -2.03
// id 45 year 2006 didn't converge with multiple guesses of starting values.
// id 77 year 2001, 2002, 2006 converged with starting value 0.44, 3.27, 0.78, 8.32

forv idval=1/111{
forv yearval=2001/2010{
use "./data/pm10data.dta", clear
qui keep if id==`idval' // city by city
qui keep if year==`yearval' // year by year
if _N>100{
scalar cff1=0.15
sum pmconc if pmconc<=float(cff1)
gen Fxnp=r(N)/_N  
gen cens=0

*****************
* CENSORED MLE *
*****************
scalar lftcff1=0.135  
scalar rtcff1=0.21
replace cens=1 if pmconc>float(lftcff1) & pmconc <=float(cff1)  
gb2lfit_mod pmconc, cdf(pmparacdf2)  censvar(cens) 

//0.31,-6.88,8.00,2.51
//
/*scalar lna0=0.31
scalar lnb0=-6.88
scalar lnp0=8.00
scalar lnq0=2.51
matrix intv = (lna0, lnb0, lnp0, lnq0)
matrix colnames intv = ln_a:_cons ln_b:_cons ln_p:_cons ln_q:_cons
gb2lfit_mod pmconc, cdf(pmparacdf2)  censvar(cens) from(intv)  */

local ahat=e(ba)
local bhat=e(bb)
local phat=e(bp)
local qhat=e(bq)
local atx=cff1
gen Fxp=ibeta(`phat',`qhat', (`atx'/`bhat')^`ahat'/(1+(`atx'/`bhat')^`ahat'))
gen propFx=(Fxnp-Fxp)/Fxnp

local cityname=city[1] 
cumul pmconc, gen(pmecdf) 
gen cutoff=pmconc*0+cff1
line pmecdf pmconc if pmconc<=0.4, sort xlabel(0(0.1)0.4) lcolor(black) graphregion(color(white)) ///
|| line pmparacdf pmconc if pmconc<=0.4, lcolor(navy) sort ///
|| line pmecdf cutoff, lcolor(black) lpattern(dot) title("`cityname': `yearval'")  ///
name(`cityname'`yearval'_CMLEfit, replace) aspect(0.8) ///
legend(order(1 "Reported Data" 2 "Estimated Counterfactual Distribution") size(4) region(c(none)) bm(zero)) ytitle("") 
graph export "./outputgraphs/app/`cityname'`yearval'_CMLEfit_alt.pdf", replace

keep in 1
if `idval'==1 & `yearval'==2003{
save "./data/PM10CMLEresults_alt.dta", replace
}
else {
append using "./data/PM10CMLEresults_alt.dta"
save "./data/PM10CMLEresults_alt.dta", replace
} 
}
}
}
